# Replication File for
# "Coalition Inclusion Probabilities: A Party-Strategic
# Measure for Predicting Policy and Politics
# Authors: Mark A. Kayser, Matthias Orlowski
# and Jochen Rehmert
# Code to replicate Table 6 and G3 and G4 (both Appendix)

# load required packages
library(stargazer);library(lmtest);library(sandwich) 

# set your directory to the folder with
# "data_eps.RDS" data file 
setwd('...')

# read data
dat <- readRDS("data_eps.RDS")

# coding of control variables 
dat$kyoto <- 0
dat$kyoto[dat$country == "austria" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "czechia" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "denmark" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "estonia" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "finland" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "germany" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "greece" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "ireland" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "italy" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "japan" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "netherlands" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "new zeland" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "norway" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "poland" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "portugal" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "slovakia" & dat$year >= 1999] <- 1
dat$kyoto[dat$country == "slovenia" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "spain" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "sweden" & dat$year >= 1998] <- 1

dat$fukushima <- 0
dat$fukushima[dat$year >= 2011] <- 1

dat$decades <- NA
dat$decades[dat$year >= 1990 & dat$year <= 1994] <- 1
dat$decades[dat$year >= 1995 & dat$year <= 1999] <- 2
dat$decades[dat$year >= 2000 & dat$year <= 2004] <- 3
dat$decades[dat$year >= 2005 & dat$year <= 2009] <- 4
dat$decades[dat$year >= 2010 & dat$year <= 2014] <- 5



dat$country_id <- as.numeric(as.factor(dat$country))

# run base-line model first to acquire row index of constant sample
summary(mod.1.cip <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                         eco_cip_gross + as.factor(decades) , 
                     data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.1.cip.se <- coeftest(mod.1.cip, vcov=vcovHC(mod.1.cip, "HC1"))[,2]


summary(mod.1.poll <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                       I(eco_polls/100)  + as.factor(decades) , 
                     data = dat[!is.na(dat$eco_cip_gross) & dat$cee == 0 ,]))
mod.1.poll.se <- coeftest(mod.1.poll, vcov=vcovHC(mod.1.poll, "HC1"))[,2]

summary(mod.1.seat <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                          I(eco_seatsh/100)  + as.factor(decades), 
                     data = dat[!is.na(dat$eco_cip_gross) & dat$cee == 0 ,]))
mod.1.seat.se <- coeftest(mod.1.seat, vcov=vcovHC(mod.1.seat, "HC1"))[,2]

summary(mod.1.all1 <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                          eco_cip_gross + I(eco_polls/100) + as.factor(decades), 
                        data = dat[!is.na(dat$eco_cip_gross) & dat$cee == 0 ,]))
mod.1.all1.se <- coeftest(mod.1.all1, vcov=vcovHC(mod.1.all1, "HC1"))[,2]

summary(mod.1.cip.environ <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                                 eco_cip_gross  + environ_worry_interpolated + as.factor(decades) , 
                       data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.1.cip.environ.se <- coeftest(mod.1.cip.environ, vcov=vcovHC(mod.1.cip.environ, "HC1"))[,2]


summary(mod.1.all2 <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                          eco_cip_gross + I(eco_seatsh/100) + I(eco_polls/100)  + as.factor(decades), 
                       data = dat[!is.na(dat$eco_cip_gross) & dat$cee == 0 ,]))
mod.1.all2.se <- coeftest(mod.1.all2, vcov=vcovHC(mod.1.all2, "HC1"))[,2]



summary(mod.1.int <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                         cab_minority * eco_cip_gross  + as.factor(decades), 
                       data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.1.int.se <- coeftest(mod.1.int, vcov=vcovHC(mod.1.int, "HC1"))[,2]


# ------------------- #
# ----- Table 6 ----- #
# ------------------- #
stargazer(list(mod.1.cip, mod.1.poll, mod.1.seat, mod.1.all1, mod.1.cip.environ, mod.1.all2),
          se = list(mod.1.cip.se, mod.1.poll.se, mod.1.seat.se, mod.1.all1.se, mod.1.cip.environ.se, mod.1.all2.se),
          omit = c("country_id", "decades"), covariate.labels = c("Minority Cabinet", "Kyoto Protocol", 
                                                                  "Quarterly GDP Growth (yearly mean)",
                                                                  "Cabinet's Mean Environmental Protection",
                                                                  "Total Greenhouse Gas Emissions/Capita",
                                                                  "Green Party in Government",
                                                                  "Green Party's Environmental Protection",
                                                                  "Green Party's gross CIP (yearly mean)",
                                                                  "Green Party's Polls", "Green Party's Seatshare",
                                                                  "Share of Population Worried About Environment"),
          add.lines = list(c("No. of Countries", length(unique(mod.1.cip[["model"]]$`as.factor(country_id)`)),
                             length(unique(mod.1.poll[["model"]]$`as.factor(country_id)`)),
                             length(unique(mod.1.seat[["model"]]$`as.factor(country_id)`)),
                             length(unique(mod.1.all1[["model"]]$`as.factor(country_id)`)),
                             length(unique(mod.1.cip.environ[["model"]]$`as.factor(country_id)`)),
                             length(unique(mod.1.all2[["model"]]$`as.factor(country_id)`))),
                           c("Country Fixed-Effects", "YES", "YES", "YES", "YES","YES" ,"YES"),
                           c("Period Fixed-Effects", "YES", "YES", "YES", "YES","YES" ,"YES") ),
          out = "table6.tex")

# ----------------------------- #
# ----- Data for Table G3 ----- #
# ----------------------------- #

# vif
library(car)
mean(vif(mod.1.cip)[,3])
mean(vif(mod.1.poll)[,3])
mean(vif(mod.1.seat)[,3])
 
# heteroskedasticity
library(lmtest)
bptest(mod.1.cip)
bptest(mod.1.poll)
bptest(mod.1.seat)
 
# outliers/influential observations
mean(cooks.distance(mod.1.cip))
mean(cooks.distance(mod.1.poll))
mean(cooks.distance(mod.1.seat))
 
# RMSE
library(Metrics)
rmse(mod.1.cip$model$eps, mod.1.cip$fitted.values)
rmse(mod.1.poll$model$eps, mod.1.poll$fitted.values)
rmse(mod.1.seat$model$eps, mod.1.seat$fitted.values)


##################################### 
#         robustness checks         #
#####################################

# run base-line model first to acquire row index of constant sample
summary(mod.1.cip <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                         eco_cip_gross + as.factor(decades) , 
                       data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.1.cip.se <- coeftest(mod.1.cip, vcov=vcovHC(mod.1.cip, "HC1"))[,2]

summary(mod.2.cip <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                         eco_cip_gross + as.factor(decades) + fukushima , 
                       data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.2.cip.se <- coeftest(mod.2.cip, vcov=vcovHC(mod.2.cip, "HC1"))[,2]

summary(mod.3.cip <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                         eco_cip_gross + as.factor(decades) + pm_cip_eco , 
                       data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.3.cip.se <- coeftest(mod.3.cip, vcov=vcovHC(mod.3.cip, "HC1"))[,2]


summary(mod.4.cip <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                         eco_cip_gross + as.factor(decades) + pm_polls , 
                       data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.4.cip.se <- coeftest(mod.4.cip, vcov=vcovHC(mod.4.cip, "HC1"))[,2]

summary(mod.5.cip <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                         eco_cip_gross + as.factor(decades) + environ_pm , 
                       data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.5.cip.se <- coeftest(mod.5.cip, vcov=vcovHC(mod.5.cip, "HC1"))[,2]

summary(mod.6.cip <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                         eco_cip_gross + as.factor(decades) + time_election , 
                       data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.6.cip.se <- coeftest(mod.6.cip, vcov=vcovHC(mod.6.cip, "HC1"))[,2]

summary(mod.7.cip <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                         eco_cip_gross + as.factor(decades) + pm_fam , 
                       data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.7.cip.se <- coeftest(mod.7.cip, vcov=vcovHC(mod.7.cip, "HC1"))[,2]


# ------------------- #
# ----- Table G4 ---- #
# ------------------- #

stargazer(list(mod.2.cip, mod.3.cip, mod.4.cip, mod.5.cip, mod.6.cip, mod.7.cip),
          se = list(mod.2.cip.se, mod.3.cip.se, mod.4.cip.se, mod.5.cip.se, mod.6.cip.se, mod.7.cip.se),
          omit = c("country_id", "decades"), covariate.labels = c("Minority Cabinet", "Kyoto Protocol", 
                                                                  "Quarterly GDP Growth (yearly mean)",
                                                                  "Cabinet's Mean Environmental Protection",
                                                                  "Total Greenhouse Gas Emissions/Capita",
                                                                  "Green Party in Government",
                                                                  "Green Party's Environmental Protection",
                                                                  "Green Party's gross CIP (yearly mean)",
                                                                  "Post-Fukushima", "PM's CIP excl. Green Parties",
                                                                  "PM Polls","PM Environmental Protection",
                                                                  "Time to Election (Years)","PM is Christian Dem.",
                                                                  "PM is Conservative","PM is Liberal","PM is Social Dem."),
          add.lines = list(c("No. of Countries", length(unique(mod.2.cip[["model"]]$`as.factor(country_id)`)),
                             length(unique(mod.3.cip[["model"]]$`as.factor(country_id)`)),
                             length(unique(mod.4.cip[["model"]]$`as.factor(country_id)`)),
                             length(unique(mod.5.cip[["model"]]$`as.factor(country_id)`)),
                             length(unique(mod.6.cip[["model"]]$`as.factor(country_id)`)),
                             length(unique(mod.7.cip[["model"]]$`as.factor(country_id)`))),
                           c("Country Fixed-Effects", "YES", "YES", "YES", "YES","YES" ,"YES"),
                           c("Period Fixed-Effects", "YES", "YES", "YES", "YES","YES" ,"YES") ),
          out = "tableG4.tex")

